////////////////////////////////////////////////////////////////////////////////////////////////////////////
//                                                                                                        //
//  -----------------------        Ebrahim Foulaadvand, 03 June 2012   -------------------------          //
//                                                                                                        //
// The programme "2dProjectile" computes the 2D trajectory of a projectile as well as its flight time,    //
// range and the maximum height in the presence of air drag (quandratic in v) in a constant gravitational //
// field g as a function of the firing angle teta0. Initial velocity v0 is fixed and given by user.       //
// Both the Euler and the Euler-Richardson algorithms are used. When you choose the algorithm,            //                                                                          //                                                                                                                 //
// Remember to inactive the section corresponding to the other algorithm.                                //
//                                                                                                        //                                                                                                         //
////////////////////////////////////////////////////////////////////////////////////////////////////////////


#include <iostream>
#include <fstream>
#include <conio>
#include <math.h>
#include <numeric>
#include <values>
#include <vector>
#include <list>
#include <algorithm>
#include <numeric>
#include<iomanip.h>

using namespace std;

double ax (double &vx,double &vy); // acceleration function for the x component
double ay (double &vx,double &vy); // acceleration function for the y component

double  tau=0.01,g=9.8,acx,acy,m=1.,c2=0.1,teta0,delteta0=M_PI/50.,v0=30,h=10.,x0=0.,y0=h,v,R,Rmax=0.,
        H,Hmax=0.,T,Tmax,tetamaxR=0.,tetamaxH=0.,tetamaxT=0.,Y,tmid=0.,ymid=0.,xmid=0.,vxmid=0.,vymid=0.;

// tau=timestep, m=projectile mass, acx(y)=instantaneous x(y) component of the projectile acceleration,
// teta0=firing angle, v0=initial velocity magnitude, h=initial height above the ground, Tmax=maximum
// flight time in the loop of teta0, Rmax=maximum range in the loop of teta0, maximum of maximum height
// in the loop of teta0.

int N=50000,t,S=25;

main()
{

vector<double> vx(N+1,0),vy(N+1,0),x(N+1,0),y(N+1,0);

// Arrays "x" and "vx" store the projectile x component of position and velocity. Similar arrays "y"
// and "vy" do the job for the y components.

ofstream trajectory("trajectory c2=0.1.plt"); // output file for the projectile trajectory.
ofstream fileR("R c2=0.1.plt"); // output file for the projectile range R.
ofstream fileH("H c2=0.1.plt"); // output file for the projectile maximum height H.
ofstream fileT("T c2=0.1.plt"); // output file for the projectile total flight time T.

x[0]=x0;     y[0]=y0;

for(int s=0; s<=S; s++){ // loop over firing angle teta0

teta0=s*delteta0;  //cout<<"teta0= "<<teta0*57.3<<endl; //getch();

H=0.; T=0.; R=0.;

// ----------- initial conditions --------------------


vx[0]=v0*cos(teta0);   vy[0]=v0*sin(teta0); //cout<<"vx[0]= "<<vx[0]<<endl;getch();

Y=20.; t=0;

/*
//------------------------ Euler-Richardson Method------------------------


while (Y>0.01){

tmid=t*tau + 0.5*tau;
ymid=y[t] + 0.5*vy[t]*tau;
xmid=x[t] + 0.5*vx[t]*tau;
vymid=vy[t] + 0.5*tau*ay(vx[t],vy[t]);
vxmid=vx[t] + 0.5*tau*ax(vx[t],vy[t]);


vy[t+1]=vy[t] + tau*ay(vxmid,vymid);
vx[t+1]=vx[t] + tau*ax(vxmid,vymid);
y[t+1]=y[t] + tau*vymid;
x[t+1]=x[t] + tau*vxmid;

Y=y[t+1];

if(y[t+1]>H) H=y[t+1];


if(teta0==10*M_PI/100.){ if (y[t]>=0) trajectory<<x[t]<<" "<<y[t]<<endl; }

t+=1;

} // end while

// The output H gives the maximum height.

T=t*tau; // T is the computed total flight time.

R=x[t]; // R ia the computed range.

*/


// ------------------------ Euler algorithm -------------------------

while (Y>0.01){

x[t+1]=x[t] + tau*vx[t];
y[t+1]=y[t] + tau*vy[t];

vy[t+1]=vy[t] + tau*ay(vx[t],vy[t]);
vx[t+1]=vx[t] + tau*ax(vx[t],vy[t]);


Y=y[t+1];

if(y[t+1]>H) H=y[t+1];


if(teta0==10*M_PI/100.){ if (y[t]>=0) trajectory<<x[t]<<" "<<y[t]<<endl; }

t+=1;

} // end while

// The output H gives the maximum height.

T=t*tau; // T is the computed total flight time.

R=x[t]; // R ia the computed range.


if(R>Rmax) { tetamaxR=teta0; Rmax=R; }

if(H>Hmax) { tetamaxH=teta0; Hmax=H; }

if(T>Tmax) { tetamaxT=teta0; Tmax=T; }


fileR<<teta0*57.3<<"  "<<R<<endl;

fileH<<teta0*57.3<<"  "<<H<<endl;

fileT<<teta0*57.3<<"  "<<T<<endl;


//cout<<"T= "<<T<<endl;//getch();

//cout<<"R= "<<R<<endl;//getch();

//cout<<"H= "<<H<<endl;getch();


}// end s

//cout<<"Hmax= "<<Hmax<<endl;

//cout<<"Rmax= "<<Rmax<<endl;

//cout<<"Tmax= "<<Tmax<<endl;

cout<<"end."<<endl;getch();

return 0;

}

//----------------------------------------------------------------

double ax (double &vx,double &vy) {

v=sqrt( vx*vx + vy*vy);

return acx=-(c2/m)*v*vx;

}

//----------------------------------------------------------------

double ay (double &vx,double &vy) {

v=sqrt( vx*vx + vy*vy);

return acy= -m*g-(c2/m)*v*vy;

}